SCANDALOUSLY PARALLELIZABLE MESH GENERATION 



D. M. BORTZ* AND A. J. CHRISTLIEB 1 " 

Abstract. We propose a novel approach which employs random sampling to generate an accu- 
rate non-uniform mesh for numerically solving Partial Differential Equation Boundary Value Prob- 
lems (PDE-BVP's). From a uniform probability distribution tt over a ID domain, we sample M 
discretizations of size N where M N. The statistical moments of the solutions to a given BVP 
on each of the M ultra-sparse meshes provide insight into identifying highly accurate non-uniform 
meshes. Essentially, we use the pointwise mean and variance of the coarse-grid solutions to con- 
struct a mapping Q(x) from uniformly to non-uniformly spaced mesh-points. The error convergence 
properties of the approximate solution to the PDE-BVP on the non-uniform mesh are superior to a 
uniform mesh for a certain class of BVP's. In particular, the method works well for BVP's with lo- 
cally non-smooth solutions. We present a framework for studying the sampled sparse-mesh solutions 
and provide numerical evidence for the utility of this approach as applied to a set of example BVP's. 
We conclude with a discussion of how the near-perfect paralellizability of our approach suggests 
that these strategies have the potential for highly efficient utilization of massively parallel multi-core 
technologies such as General Purpose Graphics Processing Units (GPGPU's). We believe that the 
proposed algorithm is beyond embarrassingly parallel; implementing it on anything but a massively 
multi-core architecture would be scandalous. 
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1. Introduction. A wide range of linear and non-linear boundary value prob- 
lems, such as the Poisson, Helmholtz, and non-linear Poisson-Boltzmann equations, 
arise in protein folding, molecular dynamics, neutral and non-neutral plasmas, acous- 
tics and antenna design, etc. In many of these applications, sparse mesh-points are 
sufficient over vast regions of the domain, while large, iterative, and parallel com- 
putations using a high density of points are necessary for obtaining accuracy in the 
vicinity of complex objects. In this paper, we propose a novel approach for generating 
mesh-points to be used in simulating solutions to Boundary Value Problems (BVP). 
The class of methods is designed to discover accurate non-uniform discretizations via 
a sparse stochastic approximation. 

State of the Art numerical methods typically have some mechanism for refining 
the underlying mesh. At some level, all adaptive methods are attempting to identify 
the optimal mesh, providing uniform accuracy over a given domain. Adaptive methods 
have been heavily employed in problems with complex geometry, where sharp geomet- 
ric features can make the simulation difficult to resolve, or in simulations where shocks 
or discontinuities can arise over time. 

Ideally, we would like an unstructured mesh that is designed to use coarse res- 
olution where the solution changes little, and fine resolution where transitions take 
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place, such as near a boundary. This is not a new idea, as this is the guiding prin- 
ciple in deciding where to apply Adaptive Mesh Refinement (AMR) for Hyperbolic 
conservation laws. Recent work in Hyperbolic conservation laws has gone as far as 
developing metrics to guide refinement which are based on an estimated error between 
the numerical solution and the weak solution |14[I15). Like with hyperbolic problems, 
the idea of working on the smallest mesh possible while maintaining a uniform ac- 
curacy is not a new topic when it comes to the area of numerical boundary value 
problems. For boundary value problems with Lipschitz sources there are a range of 
methods which can be employed to increase accuracy at boundaries where there are 
transitions. One of the simplest is to make use of Chebyshev points, which pack 
resolution near boundaries and still allow one to construct high order approximate 
inverse operators [2]. However, in a wide range of nonlinear boundary value problems, 
there are interesting issues which arise, including the formation of interior layers or 
singularities at unknown locations. In ID, if one knew where a transition layer or sin- 
gularity is located, one could design Chebyshev points to increase resolution nearby. 
Unfortunately, it may be non-trivial to identify the transition layer or develop an 
efficient implementation in higher dimensions. 

Quasi Monte Carlo and low discrepancy sequences offer an alternative approach 
to sparse mesh theory (4[ [19l [20]. A Monte Carlo method is a statistical approach to 
compute the volume of a complex domain in high dimensions. Briefly, given a complex 
domain, D Cl one encompasses the domain in a simple region D and then points are 
sampled from D and tested to see whether they lie inside of D c . The fraction of the 
points from D that lie inside D c are used to estimate the volume of domain D c . As 
the number of samples, N, go to infinity, it is well known that this method converges 
at a rate of 0(s/~N). A quasi- Monte Carlo method is formulated in the same way, 
except that instead of using randomly generated points, a pre-determined set of points, 
called a low-discrepancy sequence, is used in the estimation of the volume D c [12J. 
The Quasi Monte Carlo method converges at a rate of O(N). If one views a PDE 
from a numerical Green's function formulation, low-discrepancy sequences are one 
approach to numerically computing a solution in complex domains. These methods 
have been used in the solution of integral equations [6, 7J, and more recently, in the 
solution of PDE's \7\. The basic premise is that a good estimate of the solution may 
be obtained with relatively few samples, i.e., a sparse mesh. 

Our proposed strategy differs from all the above frameworks in that it uses sparse 
stochastic sampling to construct a non-uniform mesh-generation function Q. The 
approach is ideal for identifying the placement of mesh points to maintain uniform 
accuracy across the domain and is well suited to both finite difference and finite el- 
ement formulations of a problem. Details of the approach and associated theoretical 
framework are described in Section[2]and include a preliminary investigation of conver- 
gence and well-posedness issues. Numerical examples are then presented in Section [31 
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Section|3]contains a summary of our results as well as a discussion of future directions. 
Lastly, since this work lies at the intersection of numerical PDE's and statistics, some 
notation may be unfamiliar to some readers and accordingly, we include a glossary of 
notation in the Appendix. 

2. Approach and Theoretical Framework. Non- uniform discretizations can 
offer superior solution accuracy and convergence properties. The challenge, however, 
lies in efficient identification of a mesh which provides results superior to that of 
uniform spacing. In this section, we offer a brief overview of our proposed algorithm 
as well as the establishment of a preliminary theoretical framework. 

Consider a two-point Boundary Value Problem (BVP), e.g., the Poisson equation 

(2.1) u"(x) = f(x) s.t. u(0) = A ; u(l) = B , 

for A,BeW and x G I = (0, 1). Recall that a uniform discretization for the second 
derivative yields second-order convergence. We denote I as the domain [0,1], a dis- 
cretization of the interior I as x„ = (xi, . . . , x n ) T , and the corresponding approximate 
solution at x„ as u„ = (ui, . . . ,u n ) T . Note that frequently in the numerical analysis 
literature, the variables for spatial discretization and numerical solution are denoted 
by a superscript n. Here n is a subscript, which is consistent with notation for random 
variables and will be used in this paper. 

Next, consider a monotonically non-decreasing function Q : I — > I which is ab- 
solutely continuous on a finite number of compact subsets of I and restricted at the 
endpoints to Q(0) — 0, Q(l) = 1. The purpose of the function Q is to map the 
uniformly spaced mesh to a non-uniformly spaced one. We reserve Q° to denote the 
identity mapping Q (x) = x. 

The goal is to develop a strategy for identifying a Q such that, e.g., the approxi- 
mate solution to the Poisson problem 

u"{Q(xj) = f(Q(x)) s.t. u(Q(0)) = A ; u(Q(l)) = B , 

has convergence properties (in n) superior to a uniform spacing. We use a superscript 
to denote a uniform spacing or a variable derived from uniform spacing, e.g., x° = 

(^TT' • • • ' ^Tl) T ' u ° = ( u i> u 2> ■ ■ ■ >w°) T > Q°( x n) = x «, etc - The core of our approach 
is to identify Q via a sparse stochastic approximation. We repeatedly sample from a 
distribution P and then use pointwise statistical moments of the coarse solutions to 
generate the desired non- uniform mesh function Q. D . Naturally, different classes of 
problems call for different strategies for generating Q. Our results, however, suggest 
that a more generalizable strategy may exist. 

The rest of this section is devoted to establishing notation, framework, and im- 
plementation details. We describe notation in Section 12.11 and recall the relevant 
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development of non-uniform finite difference (FD) approximations in Section 12.21 In 
Section 12. 3[ we address how conventional consistency and stability result can be cast 
in our framework. In Section [2.41 we motivate the construction of Q as well as present 
the algorithms for the test cases in Section ([3j) . Lastly, Section 12.51 is devoted to a 
description of the sampling strategy. 

2.1. Notation. Let X n = (Ai, A2, ■ ■ ■ , A n ) be a random vector where the ele- 
ments are independent and identically distributed (i.i.d.) with A& ~ P for a given 
probability distribution P. The points in the realizations must be sorted before 
being used to construct an approximate solution. It is statistical convention to 
denote a sorted sequence using subscripts with parenthesis. Accordingly X(„) = 
X(2), . . . , X(„)) T denotes a sorted random vector with elements from X„ where 
x (k) ~ (™)-P j (l - P) k ~ 3 - The derivation for this distribution comes from order 

statistics and we direct the interested reader to 

Definition 2.1. Let U : I" — > W l be a mapping taking a monotonically increasing 
sequence of n distinct points in I and solving a given BVP on those points. Note that 
U represents the action of numerically solving a PDE and in some of the examples 
below, we discretize the second derivative using a simple three-point stencil. In several 
cases, we implemented higher order stencils (results not presented), and while they 
do provide better convergence results, the conclusions of this work do not change. 

Definition 2.2. Let p be a function taking a point £ £ I and a random vector of 
length n, and mapping them to a single random variable 

(2.2) p(^X {n) (P))=E K [{U(X {n) (P))} K=k \X {k) =£] . 

The operator E#- denotes expectation with respect to a uniform distribution on 
{1, . . . , n} where the distribution of the index random variable K and {-}k denotes 
the Ath element of a vector. We note that this function returns a random variable 
for each £. 

Definition 2.3. Let the pointwise mean of p be defined for £ € I as 

(2.3) M (0 = Ep [Ex [{UpL W (P))} K=k \X [k) =(]}. 

To appreciate the utility of this definition for /i, consider that this allows pointwise 
computation of the expected solution value for an arbitrary point £ € I by conditioning 
the expectation on £ being one of the points in a realized mesh. The expectation 
on an index K was constructed after much consideration and it is productive to 
reflect on alternative choices. An alternative definition of p could include information 
(via an interpolation) from solutions where £ is not one of the mesh points. We 
found, however, that this reduced the efficiency of our algorithm both theoretically 
and practically. In addition to increasing the complexity of the analysis, adding an 
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interpolation step to the consistency and stability results (Section 12. 3[) induced a 
reduction in the speed of solution convergence. 

Definition 2.4. Let the pointwise variance of p be defined for £ G I as 

(2.4) W (0 = V P [E K [{C/(X ( „)(P))} K=fc |X (fe) = £]] , 

where Vp denotes variance with respect to P, E^ denotes expectation with respect to 
14{1, . . . ,n}, the distribution of the index random variable K, and {-}k denotes the 
Kth element of a vector. We also let the average variance over I be defined as 

(2.5) v = p/ v ^ [{^(x ( „)(p))} K=fc |A: (fc) = e]]de. 

In the examples presented in Section [3l we choose P to be a uniform distribution 
on I. There is nothing in the following development, however, that would fundamen- 
tally change for an arbitrary P. 



2.2. Conventional Non-uniform Finite Difference. In the example prob- 
lems of Section [31 we will employ finite differences. There are first and second deriva- 
tives in these problems, which must be discretized on a nonuniform mesh. For the first 
derivative, we select the simple (though unstable) centered difference discretization. 
For the second derivative, however, we create a non-uniform centered difference ap- 
proximation to the differential operator. For uniform mesh x° and any non-uniform 
mesh function Q, we can Taylor expand the solution u at mesh-points k + 1 and k — 1 
around node k 

u k +i =u k + h k+1 u' k + hil^u'l + 0{hl_ x ) , 
u k -i = u k - h k u' k + -h%u'l - 0{h\) , 

where h k — Q(x k ) — Q(x k _ 1 ) is the step-size and u k is the solution evaluated at Q(x k ). 
This leads to a Method of Undetermined Coefficients problem 

a k Uk+i + b k u k + c k u k _i = (a k + b k + c k )u k + {a k h k+l - c k h k )u' k 

+ \( a kh 2 k+1 +c k h 2 k )u'l 1 

with the algebraic equations 

a k + b k + c k = 0, a k h k+ i ~ c k h k = 0, -{a k h\ +1 + c k h k ) = l. 
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The values of a^, b^, and c^, yield the discretized version of (|2.1[) . ^4q( x o )U„ = /q( x o j 
where 4q( x o) is a tridiagonal matrix 



(2.6) 
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/q(x°) is the forcing function is evaluated at the n mesh-points, and u n is the ap- 
proximate solution at the Q(x°) mesh-points. Relating this to the notation above, 
we have that u„ = L/(Q(x°)) = ^q( x o)/q( x °)- 

We also note that in what follows, we consider a randomly sampled and sorted 
mesh X(„)(P) and the corresponding approximate differential operator is denoted as 



A* 



)(P) 



2.3. Consistency and Stability. For a conventional finite difference discretiza- 
tion, we would consider the error E in the solution 



\E(Q,<) 



U (Q(x°))-L/(Q(x°))|| 
0) (A Q(x o )W (Q(x°)) 



< 



A' 1 

A- 1 



no 



\ T Q«)\ 



which is bounded above by the spectral radius of the inverse of the finite difference 
operator ^4q( x o) an d a truncation error tq( x o). For the non-uniform three-point- 
stencil approximating the second derivative, the truncation error is ©(max^ \hk\)- 

For our development, however, we consider a probabilistic version of this error, 
with the following conditions. 

CONDITION CI. For a given P, the spectrum of A^ is bounded in [0, 1]. 

Condition C2. For a given P, the truncation error induced by a finite difference 
approximation to the second derivative is first order in the largest step-size h. 

Moreover, we go so far as to propose that Conditions CI and C2 are in fact true 
for all classes of BVP's. In support of Condition CI, the red circles in Figure [27X1 
depict the maximum and minimum eigenvalue magnitudes for m = 12000 realizations 
of X( n ) (P) for P a uniform distribution on I. For comparison, we have included the 
eigenvalue bounds for a uniform partition (blue squares). In support of Condition C2, 
the red circles in Figure [2~2l depict (for each n) the maximum step-size in 1000 meshes 
sampled from P, a uniform distribution on I. This illustrates that as n increases, the 
maximum step-size is decreasing as well, albeit slowly in comparison with step-size in 
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Figure 2.1. Red circles are the magnitude of upper and lower bounds for eigenvalues of 
j(p)- F° r each n, m = 12000 random grids were sampled from a uniform distribution and 

eigenvalues of the corresponding A~^~ were computed. For comparison, the blue squares are the 

magnitudes of the upper and lower bounds for the eigenvalues of A~\ . 




Figure 2.2. Red circles illustrate the convergence of the largest step-size in randomly generated 
meshes of length n. For each n, m = 1000 meshes X( n j(P) were generated for P a uniform distri- 
bution on I. The vertical axis is the maximum step-size between and two nearby points in X( n )(P). 
For comparison, we also include the step-sizes for a uniform partition of the domain (blue squares). 



a uniform partition (blue squares). 

Theorem 2.5. Under Condition CI and C2, the expected error converges point- 
wise to zero. 

Proof. Let be a vector of zeros with a one in the fcth element and X^v^fc 
denote the vector (-X"(i), ■ • ■ , Xt^-i) , £, Xr^%\, . . . , Xr n \) where the fcth element of X(„) 
is replaced by £. Consider the pointwise error between the the analytical solution u 
and mean of the coarse-mesh solutions fi 

IKO - M£ U(X {n) (P)))\\ = \\u(0 - E P [E K [{U(X (n) (P))} K \K : X {K) = £]] || 

= ||E P [E K [u(0 - {U(X (n) (P))} K \K : X (K) = £]] || 

= ||E P [E K [u(Q - el(U(X (n) (P)) \K : X (K) = £]] || 

< E P [E K [||e£(u(X ( „ UiJS -) - U(X (n) (P))\\ \K : X {K) = £]] 
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We note that A^.,(p)u(X( n ),£,k) — fx (n) is the truncation error induced from a 3- 
point non-uniform centered-difference stencil and converge as ©(max^i,. ,, )n _i \Xu\ — 
(Figure also depicts this behavior). Therefore, assuming Conditions CI 
and C2, we can conclude that the error converges pointwise. □ 

2.4. Criteria Motivation. In Section [3J we consider problems containing sec- 
ond derivatives as well as nonlinear functions of first derivatives. We create Q de- 
pending upon the specific problem class under consideration. In one case, Q is based 
upon the derivative of the mean of the solutions fi, while in another we consider the 
product of the variance and the the second derivative of /i. In all cases, however, Q 
is created using the statistical moments of the sampled sparse-mesh solutions and the 
results suggest that a more generalizable strategy may exist. 

For the problems with second derivatives (equation (|3.1[) in Subsection 13.11 and 
equation (|3 . 6[) in Subsection 13. 3p we define Q as 



Q(x) 



?l(0 



(x). 



qi(x)= [ y/\li'(t; E/(X (n) (P))|d£, 



where 



and the superscript —1 is an inverse function operator. Essentially, this definition 
will pile up points in regions with a steep solution in an effort to provide higher order 
accuracy for the nonuniform second derivative discretization. 

For the problem with a third power of the first derivative (equation (|3.3p in 
Subsection l3.2[) . we define Q as 



where 



Q(x) 



92(0 



92(1) 



(x). 



qi{x) = / ^"(eU(X (n) (P)) 2 v(eU(X (n) (P)) 3 d^ 



and v is from Definition (|2.4j) . 

To motivate these mappings, consider a least squares cost function J : M — > 
R which compares (pointwise) a proposed solution to the actual solution u at an 
arbitrary point £ € I. We then consider the random variable J(/x(£, X( n )(P)) and 
Taylor series expand the expected value of J around a local minimum for some 



where 



and 



E P [j(/x(£,X ( „)(P))] a E P [J(/x*(0) + -A + J'- 



2\ 



Jx = J'(M*(0)(MC,x ( „)( p )) - 



j 2 = j'V(0)(M(e,x (n) (i 3 ))-M*(0) • 

Note that since /it* is assumed to be a local minimum for J , Ji = and that in J2, 
J"(/z*(£)) = 1. The expected cost function can thus be approximated by 

(2.7) E P [J(/i(£,X (n) (P)))] «Ep[J( M *(£))]+Ep [(^X;^?))-^^)) 2 ' 

According to Theorem 12. 51 lhrin^ooEp [J(a**(£))1 = pointwise. Most importantly, 
the second term in ()2.7p is the variance of /x(£, X(„)(P)). The moments of the sparse 
sampled solutions, therefore, can offer insight into the actual solution. 

We hypothesize that the reason q\ (x) works well is that the // may converge faster 
than fi. A full investigation will require establishing the appropriate Sobolev space, 
though it is not immediately clear how to proceed. 

We also hypothesize that the function qi (x) works well because the second deriva- 
tive (when cast as the local curvature) is inversely proportional to the local variance 
of a random variable (a result which is well known in the semi-parametric nonlinear 
regression literature [23]). Essentially, while the fi" may not converge quickly, the 
product fi"v does. We also found that multiplication by an extra v dramatically im- 
proves the computed Q, though an explanation is not immediately clear. A deeper 
understanding of the spectrum of A X(n) (p) and how it depends upon the choice of P 
will be essential to explaining the efficiency of q% (x). We plan to explore both of these 
issues in a future paper. 

2.5. Sampling Strategy. The proposed methodology relies on random sam- 
pling and sorting a sampled vector from a distribution. For standard sampling 
problems, the confidence intervals for variance estimation converge as 0(y/m) for 
m samples. In the absence of a full theoretical analysis, we offer numerical evidence 
supporting our accurate computation of the variance. Note that in all the following 
simulations, we use the first example problem presented in Section [3l the singular 
two-point BVP ([3~Tj) - (j3T2]) . 

The relationship between the mesh size n and the number of samples m is non- 
trivial, and Figure 12.31 illustrates this by depicting the error in v (relative to v com- 
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v(U(T)) error vs. m and n 




Figure 2.3. log ln of the error in the computation of v (sampling from a uniform distribution 
on I) as a function of m and n . Note the general downward trend along both the m and n axes. 

Samples needed for 3 digits of v(li(l)) accuracy 
1500^ 

o o 

■|iooo u O q 
4j O 

1 500 O 

O 

o 
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n 

Figure 2.4. For each n, the vertical axis reflects the number of samples needed to compute 
the variance with 3 digits of accuracy relative to v (sampling from uniform distribution onl) with 
m = 3000. 



puted with m — 3000 sampled from a uniform distribution on I) for a range of n and 
m values. For a given n, though, we do note that the error in the v computation is 
decreasing. In Figure [2"^tl we depict the number of samples of vector size n which are 
needed to ensure three digits of accuracy in estimating the variance. Since the number 
was consistently below 1000 over a range of n, we let m = 15000 in all subsequent 
simulations (unless otherwise specified). 

We also note that figures with similar conclusions can be generated for the other 
problems described in Section [3J 

3. Numerical Simulations. In this section, we consider three classes of BVP's 
and illustrate the effectiveness of our approach in identifying a non-uniform mesh 
spacing which yields a solution accuracy superior to that of uniform spacing. The 
examples in Sections I3.HI3.2I were chosen because the solutions exhibited changes in 
smoothness such as boundary layers and a discontinuity in the derivative, i.e., a local 
effect. The example in Section [3.3l was chosen because the solutions exhibited varying 
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Solution to Singular BVP for a = 0.735, (5 = 10 
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Figure 3.1. Analytical solution to singular two-point Boundary Value Problem in (3. l}) - (3.21il 
for f(x, u) = 0x a +P- 2 ((a + 0-1) + 0x^)u, u(0) = 1, and u(l) = e. 

multi-scale oscillations over the domain, i.e., a global effect. 

All software development was done in Matlab and in particular, parallel for- 
loops (parfor) were used whenever possible. For all simulations, we let m = 15000 
unless otherwise specified. A copy of the software is available upon request to the 
corresponding author. 

3.1. Singular Boundary Value Problem. The following singular two-point 
Boundary Value Problem 



encompasses a large class of singular BVP's where a € (0, 1) and A, B are finite 
and constant. For (x,y) € {[0,1] x R} such that f(x,y) is continuous, df/dx ex- 
ists and is continuous and df/dy > 0. Many researchers have successfully imple- 
mented strategies to solve classes of singular BVP's using uniform finite difference 
|13j . Rayleigh-Ritz-Galerkin [9J, projection [2T], and collocation [55] schemes. 

We consider the specific problem presented in [3J for f(x, u) — (3x a+l3 ~ 2 ((a + (3 — 
1)+ /3x^)u, A = 1, B = e and with exact solution u[x) = exp(x^). Depicted in Figure 
(|3.ip is the solution for a = .85, and j3 = 10; note the strong boundary layer affect 
near x = 1. 

Inspired by the fact that a non-uniform mesh provides superior accuracy for this 
class of problems (as investigated in [31 [S]), we consider the mapping Q a {x) — x 1_Q . 
We choose a = 0.735, Qo. 735(2;) = x 1 ^ - 735 because for n = 100, this Q0.735 yields 
2.4 orders of magnitude more accuracy in the solution than Q . The dotted curve in 
Figure [3~3l depicts the corresponding plot of Q0.735 an d should thus be considered a 
gold standard, with which we will compare our computed mapping. 

For this problem, we first consider m = 15000 samples of length n = 10 and 



(3.1) 
(3.2) 



(x a u'Y = f(x,u); < x < 1 
tt(0) = A; u(l) = B 
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Solution Density: m = 15000, n = 10 




0.5„ L 
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Figure 3.2. Contour plot of the probability density of p(x, X( n j (P)) for m = 15000 sampled 
meshes of length n = 10 from mesh mapping Q° and P a uniform distribution on I. 



solve (|3.1jl - l|3.2jl . Figure 13.21 depicts a contour plot of the probability density of 
X(„)(P)), where P is a uniform distribution on I. We note that many of the 
sparse solutions have relatively large values for x near 3/4, but we are at a loss to 
explain this phenomenon. Figure 13.31 depicts our computed Q along with the uni- 
form and analytically optimal mappings and Figure GQ] depicts the error convergence 
properties for Q and Q°. We draw attention to the fact that for mesh sizes less than 
100, our mapping results in second order convergence, even though the discretization 
is first order. Indeed, the error is at times, marginally better than Q a . We also note 
that while Q and Q a are second order accurate, our non-uniform mapping asymptot- 
ically converges with first-order accuracy. Lastly, Figure 13.5a depicts solutions using 
n = 20 mesh-points from the uniform and computed optimal mappings. The uniform 
mesh creates a solution which is consistently above the analytical solution, while the 
solution computed from Q is substantially more accurate for small numbers of mesh 
points. While the difference between the solution is visually small, it is almost an 
order of magnitude more accurate and Figure [33b depicts the solution residuals. 

3.2. Hamilton-Jacobi Equation. An important example problem which may 
develop a singularity at an internal location is the time-invariant Hamilton-Jacobi 
(HJ) equation, 

(3.3) u + H(u x )=f ;0<x<l, 

(3.4) u(0) = u(l) 

where we assume period Dirichlet boundary conditions. We consider 

(3.5) H{p) = ^ , f{x) - - sin(4* - £)) + cos 2 (4* - £)) , 



4' 
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Figure 3.3. Comparison of different mesh mappings. The solid line is uniform spacing Q°(x), 
the dotted line is the analytically optimal non-uniform spacing z( 1_c "), and dot-dashed line with the 
circles is the computed Q(x). 
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Absolute Value of Residuals for n = 10 
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Figure 3.5. a) Solution using n = 10 mesh points. The solid line is the analytical solution. 
The 'X' and 'O' symbols are the solutions generated using the uniform (Q°) and computed (Q) mesh 
point mappings, b) The absolute value of the pointwise residuals for both solutions. 
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Solution to Time-Independent Hamiltion-Jacobi 
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Figure 3.6. Analytical solution to time-independent Hamilton- J acobi equation in (3. Sty - HT^ 
for H{p) = £ and f(x) = - \em(%(x - f ))| + cos a (7r(* - f )). 

Solution Density: m=15000, n=10 
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Figure 3.7. Contour plot of the probability density of p(x, X( n ) (P )) for m = 15000 sam- 
pled meshes of length n = 10 for Hamilton- J acobi problem in i3.3\ )- i3. 51 ). where P is a uniform 
distribution on I. 



which yields a solution u(x) — — | sin (-zr (x — j)) \ with a discontinuity in the deriva- 
tive at j . The analytical solution on I is depicted in Figure 13.61 

For this problem, we first consider m — 15000 samples of length n = 10 and solve 
(|3.3p - (|3.5p . We discretized u x using centered differences and then solved this nonlinear 
PDE by minimizing the norm of the residual (u + H(u x ) — /) using a Quasi-Newton 
method with BFGS update. Figure 13.71 depicts the contour plot of the probability 
density of X( n )(P)) where P is a uniform distribution on I. 

Figure 13.91 depicts our computed optimal Q along with the uniform and analyti- 
cally optimal mappings while Figure 13.101 depicts the convergence properties for the 
different mappings. We note that computed mapping will result in many mesh-points 
near 7r/4, which is exactly the location of the derivative discontinuity in the solu- 
tion. Figure 13.101 depicts the convergence in error for uniformly and non-uniformly 
spaced points as specified by the mappings defined in Figure [XU We note that our 
strategy offers a substantial improvement in solution accuracy for small numbers of 

14 



Solution Values: m=15000, n=10 
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Figure 3.8. Point cloud of solution values for m = 15000 sampled meshes of length n = 10 
for Hamilton- Jacobi problem in 13.36 - 13.56 . 
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Figure 3.9. Mapping functions Q and Q° . Note how there will be a large percentage of the 
points at 7r/4. 

mesh-points; it even overcomes the inherent instability of centered-differences. To 
illustrate the improved accuracy for small n, Figure 13". 1 II depicts a) the solutions gen- 
erated by Q and Q° as well as b the absolute value of the pointwise residuals for both 
mesh-point mappings. 

3.3. Oscillatory Boundary Value Problems. The following problem was se- 
lected from [18] specifically because it illustrates a situation for which our approach 
does not perform well. Consider the specific Poisson problem 

(3.6) u"{x) = -20 + ^<t>"(x) cos(^(x)) - ^'{x)) 2 sm{^{x)) 

(3.7) u(0) = 1 ; u(l) = 3 

with a — 0.5. The solution with Dirichlet boundary conditions u(0) — 1, = 3 is 



u(x) = 1 + 12x - 10x 2 + asm(4>(x)) 
15 



Error Convergence For Each Mapping 
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Figure 3.10. Error convergence for uniformly and non-uniformly spaced points where the 
non-uniform spacing is determined by the mapping in Fiaure Y3.9\ 
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Figure 3.11. a) Solution to the HJ equation (3.31i l- i3.91i l with n = 40 mesh-points using the 
computed non-uniform mapping circles 'O' and uniformly spaced mesh-point 'X'. b) The absolute 
value of the pointwise residuals for both solutions. 

This analytical solution is depicted for <j>(x) = 57ra; 3 and 4>(x) — 207ra; 3 in Figures 
13.12b and 13.12b . respectively. 

As with the other problems, we first consider m — 15000 samples of length n = 10 
and solve (|3 . (|3 . 7[) . Figure [5.131 depicts contour plots of the probability density of 
p(£;X( n )(P)), where P is a uniform distribution on I, for 4>(x) = 5irx 3 and <f>{x) = 
2Qttx 3 in Figures 13.13b and 13.13b . respectively. We note the relatively large total 
variance and that higher frequency forcing oscillations appear to correlate with larger 
variance and plan to investigate this in a future article. 

Figures r3.14h ) and 13. 14b depict our computed Q for the two <fi(x) functions along 
with the uniform mapping using n — 10. Figure [3.151 depicts the convergence in error 
for uniformly and non-uniformly spaced points as specified by the mappings defined 
in Figure [3.141 Here, the Q computed by our approach is noticeably worse than using 
a uniform spacing. 
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Solution to Poisson equation with <|>(x)=5jtx 
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Figure 3.12. Analytical solution to Poisson equation from with a) <j>(x) = 57rx 3 and b) <j>(x) 
20nx 3 from JTEj . 



Solutions: m = 15000, n = 10 



Solutions: rn = 15000, n = 10 



30 
20 
= 10 

I 

3 

? -10 
-20 

-30- 



0.5 
Domain x 



0.1 

0.08 

0.06 

0.04 

0.02 





400r 



200 



-200 



-400„ L 



b) 



0.5 
Domain x 



0.1 

0.08 

0.06 

0.04 

0.02 





Figure 3.13. Contour plot of the probability density of p( x, (P )) for m = 15000 sampled 
meshes of length n = 10 for the Poisson equation defined in i3. 6V - H3. 71) , where P is a uniform 
distribution on I. 



4. Summary. 

4.1. Conclusions and Future Work. In this article, we have proposed a strat- 
egy for identifying non-uniform mappings to improve solution accuracy when using a 
smaller number of points. We find that our strategy performs well when faced with 
local phenomena, i.e., when there is a subset of the domain which would benefit from 
a higher density of mesh-points. With behaviors which are fundamentally global, such 
as oscillations, our approach fares no better than with uniformly placed points. 

The immediate goal for the next article is to apply our approach to higher dimen- 
sional systems where we envision the benefit of the 0{y/m) convergence (independent 
of n) will become even more pronounced. With the recent release of software to enable 
calling CUDA (Compute Unified Device Architecture [10J) from Matlab, we also plan 
to develop an implementation on a GPGPU. Among the open theoretical questions we 
plan to investigate in followup articles, the most pressing is how to formalize taking 
the information in the point cloud of sparse-mesh solutions and construct Q. We en- 
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Comparison of Mapping functions Comparison of Mapping functions 
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Figure 3.14. Mapping functions Q and Q° for a) <f>(x) = 5ttx 3 and b) <f>(x) = 20nx i . 




vision development of a criteria, which is based in part on the operator discretization 
strategy and part on the form of the PDE. 

Our algorithm could also greatly benefit from connections to the experimental 
design literature. There are criteria designed to aid in choosing optimal sampling 
locations, which could provide a more efficient computation of Q [1]. Our original 
vision for this work involved casting P^ 1 as the mapping Q and developing an inverse 
problem strategy to identify an optimal P* via optimizing v(P). For this paper, the 
optimization was not nearly as efficient as simply computing Q. There are, how- 
ever, stochastic approximation algorithms (such as Kiefer-Wolfowitz [T6j ) which are 
designed to optimize stochastic functions and there is significant literature within 
the control theory community on this topic |17] . This would allow reduction in the 
number of samples needed to compute the total variance. We also plan to investigate 
an improved strategy for the function approximation choices (such as how many n 
are needed). A long-term goal is to investigate how our approach could be used to 
update a given set of mesh-points to improve accuracy. Lastly, there are many oppor- 
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tunities to extend the theoretical analysis. In particular, the function theoretic basis 
for the finite element approach could prove useful in establishing convergence the the 
appropriate Sobolev space. Lastly, the rich literature of random matrix theory will of 
course improve our knowledge of the distribution of eigenvalues for the approximate 
finite difference operator. 

4.2. Discussion. The widespread availability of massively multi-core architec- 
tures such as General-Purpose Graphics Processing Units (GPGPU's) provides a 
unique opportunity for investigating non-traditional and non-conventional approaches 
to solving computational problems. It is an opportunity to completely re-imagine 
traditional numerical analysis algorithms. For example, in [8], the authors propose 
a parallel ODE solver which would be particularly slow on single-core CPU's. How- 
ever, when implemented on a GPGPU, it can realize spectacular improvements in 
computational efficiency and accuracy. Similarly, the approach presented above rests 
on the notion that future trends in architecture design will produce processors with a 
large numbers of computational cores and relatively modest memory resources. Ac- 
cordingly, we envision that the implementation of the sampling step on a GPGPU 
will yield substantial speedups in efficient mesh generation, since the algorithm will 
theoretically scale with the number of computational cores. Like the parallel ODE 
solver in [5], a single-core implementation would be completely unacceptable. We 
assert that algorithms in this class are beyond embarrassingly parallel; implementing 
it on anything but a massively multi-core architecture would be scandalous. 
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Appendix A. Notation. 

Aq^t.) Matrix generated by a three point stencil (potentially non- uniform) dis- 
cretizing the second derivative in a Poisson equation. The subscript indi- 
cates the set of points to use in the discretization. 

/q(x") Vector generated by evaluating the forcing function in a Poisson equation. 
The subscript indicates the set of points to use in the discretization. 

I Domain for the BVP. In this work, I = (0, 1) and I = [0, 1]. It is also the 

probability sample space for X. 

m Number of sampled meshes. 

ft Pointwise mean of p, which allows computation of the expected value of 

the approximate solutions at a given point £ g I, given that £ is a point 
in the mesh. 

n Number of points in a given mesh. 

p Function mapping a random vector of length n and a point £ g I to a 

random variable. 
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P Probability distribution. 

Q Mesh function which maps I — > I. 

T Q(x n ) Taylor Series truncation error from using a three point stencil to discretize 
the second derivative on a (potentially non-uniform) mesh. The subscript 
indicates the set of point to use in the discretization. 

u Analytical solution to a BVP. 

Uq n-dimensional numerical approximation to u. 

U Function which maps I" — > W 1 taking a mesh and solving a given BVP 

on those mesh points. 
v Pointwise variance of p, which allows computation of the variance of the 

approximate solutions at a given point £ £ I, given that £ is a point in the 

mesh. 

v Average of variance v over I. 

x" Mesh of n uniformly spaced points in I. 

X(P) Random variable with distribution P. 

X„ Vector of n i.i.d. random variables. X„ = (X\, X2, ■ ■ ■ , X n ) T . 

X(„) Sorted vector of n i.i.d. random variables. X(„) = (-^(i)> ^(2)1 • • • j^(n)) T - 

£ Arbitrary point in I. 
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